% This calibration assumes the benchmark cost function:
% cost(t) = c1*(Z(t))^gamma1 + c*(Rz-Rx)^gamma2

% To reatain a four-parameter calibration, set a value for
% gamma_1 from the literature. 

% Calibration procedure is robust to mutually profitable deviations:
% Having calibrated the cost function and the fee to match the targets at Rz
% = 5%, check whether Rz = 5% is indeed optimal by calculating iso-profit-fees 
% for other Rz's and checking whether a naive agent would
% prefer a deviation to those. 


% FOUR FREE PARAMETERS
% 1) c1
% 2) c2
% 3) gamma2
% 4) phi (eqm. fee expressed as a percentage of the agent's WTP fmax)

% FOUR MOMENTS
% 1) Rz = 5% maximises profits,
% 2) total cost ratio, measured as ratio of means, is equal to a specified
% target,
% 3) ratio of admin to investment costs is equal to a specified target,
% 4) markup at the optimum is equal to a specified target. 


clear
load('recalibrated_data_with_eval.mat')

%precision states, in percentage terms, how close I want to be to the
%cost-related target values
precision = 0.00001;

%fineness states how fine is the grid for c_2
fineness = 1000;


% SPECIFY TARGETS HERE:
% Rz that maximises profits
THETAtarget = 0.05;
% Total cost ratio
totalcosttarget = 0.005;
% Ratio of admin to total cost
adminratiotarget = 0.50;
% Markup at the optimum
markuptarget = 0.2;


% HERE WE GO!
totalcount = 0;
conditionscount = 0;
successcount = 0;
moveup = 0;
movedown = 0;


fileID = fopen('res_Firm_Calibration_2021.txt','wt')

rho = mydata.rho;
delta = mydata.delta;
effhh_size_ = mydata.effhh_size;
avgIncome = mydata.avgIncome;
meanhhs = mydata.meanhhs;
R = mydata.R;
zpen_eval = mydata.zpen_eval;
alpha = mydata.alpha;
survival = mydata.survival;

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_naif');
avgZ_naif_target = mydata.(suffix);
suffix = strcat('fmax_',num2str(THETAtarget*10000));
fmax_target = mydata.(suffix);

gamma1 = 0.5;

gamma2_0 = 5.60;
gamma2_1 = 5.90;

% Given gamma2, use the adminratiotarget and totalcosttarget to select an appropriate value for c2. 
c2_0 = (1-adminratiotarget)*totalcosttarget*mean(avgZ_naif_target) / (1+THETAtarget-R)^(gamma2_0);
c2_1 = (1-adminratiotarget)*totalcosttarget*mean(avgZ_naif_target) / (1+THETAtarget-R)^(gamma2_1);


% Given gamma1, use the adminratiotarget and totalcosttarget to select an appropriate value for c1.
for c1 = adminratiotarget*totalcosttarget*mean(avgZ_naif_target) / mean(avgZ_naif_target.^(gamma1))
    for gamma2 = gamma2_0:0.05:gamma2_1
        for c2 = linspace(min(c2_0,c2_1), max(c2_0,c2_1), fineness*500)
            totalcount = totalcount+1;
            
            % Check the cost conditions for Rz = 5%
            THETA = THETAtarget;
            avgZ_ = avgZ_naif_target;
            
            costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
            
            ratioatmean = mean(costs) / mean(avgZ_);
            %meanratio = mean(costs ./ avgZ_);
            adminratio = mean(c1*(avgZ_).^(gamma1)) / mean(costs);
            
            if abs((ratioatmean - totalcosttarget)/totalcosttarget) <= precision && abs((adminratio - adminratiotarget)/adminratiotarget) <= precision
                conditionscount = conditionscount + 1;
                
                % Calibrate the fee so that markup is on target. 
                phi = 0;
                eq_fee = phi*fmax_target;
                
                while (eq_fee/mean(costs) - 1) < markuptarget
                    phi = phi + 0.005;
                    eq_fee = phi*fmax_target;
                end
                
                eq_fee = min(eq_fee, fmax_target);
                
                % Calculate the resulting firm profits (PV)
                profits = ones(1,71)*eq_fee - costs;
                ProfitPV = profits(1);
                for t = 2:71
                    ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
                end
                maxprofit = ProfitPV;
                
                % Given the fee and profits at the target Rz*, calculate iso-profit-fees for other Rz's   
                for THETA = mydata.THETA_grid
                    if abs(THETA - THETAtarget) < eps
                        suffix = strcat('isofee_',num2str(THETA*10000));
                        workspace.(suffix) = eq_fee;
                        
                    else
                        suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
                        avgZ_ = mydata.(suffix);
                        
                        costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
                        
                        fee = 100;
                        profits = ones(1,71)*fee - costs;
                        ProfitPV = profits(1);
                        for t = 2:71
                            ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
                        end
                    
                        while maxprofit > ProfitPV
                            fee = fee + 20;
                            profits = ones(1,71)*fee - costs;
                            ProfitPV = profits(1);
                            for t = 2:71
                                ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
                            end
                        end
                        
                        suffix = strcat('fmax_',num2str(THETA*10000));
                        fmax = mydata.(suffix);
                        
                        suffix = strcat('isofee_',num2str(THETA*10000));
                        workspace.(suffix) = min(fee-10,fmax);
                    end
                end
                
                % Armed with iso-fees, calculate consumer's PERCEIVED
                % utility from different Rz-contracts. 
                if rho == 1
                    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
                    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
                else
                    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
                    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
                end
                
                lifetime_iso_grid = [];
                for THETA = mydata.THETA_grid
                    suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
                    avgZ_TC = mydata.(suffix);
                    suffix = strcat('avgX',num2str(THETA*10000),'_TC');
                    avgX_TC = mydata.(suffix);
                    suffix = strcat('avgTC',num2str(THETA*10000),'_TC');
                    avgTC_TC = mydata.(suffix);
                    
                    suffix = strcat('isofee_',num2str(THETA*10000));
                    isofee = workspace.(suffix);
                    avgTC_TC = avgTC_TC - ones(1,71)*isofee;

                    instant = zeros(1,71);
                    for t = 1:71
                        if rho == 1
                            instant(t) = effhh_size_(t) * log(avgTC_TC(t)/effhh_size_(t));
                        else
                            instant(t) = effhh_size_(t) * ((avgTC_TC(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
                        end
                    end
    
                    bequest = zeros(1,71);
                    for t = 1:71
                        cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
                        if rho == 1
                            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
                            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
                        else
                            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
                            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
                        end
                        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
                    end
    
                    lifetime = instant(1);
                    for t = 2:71
                        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
                    end
                    
                    lifetime_iso_grid = [lifetime_iso_grid lifetime];
                    
                    if abs(THETA - THETAtarget) < eps
                        lifetime_at_target = lifetime;
                    end
                end

                % Finally, check the equilibrium condition:
                if lifetime_at_target == max(lifetime_iso_grid)
                    successcount = successcount + 1;
                    fprintf(fileID, ['**************\r\n']);
                    fprintf(fileID, ['c1 = %f, gamma1 = %f, c2 = %24.22f, gamma2 = %f \r\n'],[c1,gamma1,c2,gamma2]);
                    fprintf(fileID, ['The cost ratio AT MEAN is %13.10f \r\n'],ratioatmean);
                    %fprintf(fileID, ['At the same time, MEAN cost ratio is %f \r\n'],meanratio);
                    fprintf(fileID, ['The ratio of admin to total costs is %f \r\n'],adminratio);
                    fprintf(fileID, ['Fee associated with requested markup is %f coming as a result of phi = %f \r\n'],[eq_fee,phi]);
                    fprintf(fileID, ['SUCCESS!!!!!! \r\n']);
                    fprintf(fileID, ['Deviation-proof, optimal contract has Rz = %f. \r\n'],THETAtarget);
                    
                else
                    [value, id_pref] = max(lifetime_iso_grid);
                    if mydata.THETA_grid(id_pref) > THETAtarget
                        moveup = moveup + 1;
                    else
                        movedown = movedown + 1;
                    end
                end
            end
        end
    end
end


fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Total runs %i \r\n'],totalcount); 
fprintf(fileID, ['Cost conditions are met in %i runs \r\n'],conditionscount); 
fprintf(fileID, ['Rz of 5 per-cent maximises profits in %i runs \r\n'],successcount);
fprintf(fileID, ['Otherwise, agent prefers a more expensive contract in %i runs \r\n'],moveup);
fprintf(fileID, ['And a cheaper one in %i runs \r\n'],movedown);
fclose(fileID);
disp(['DONE!!'])
    















